R para Ciencia de Datos en Salud:
Análisis Descriptivo e Inferencia Estadística
Percy Soto-Becerra M.D., M.Sc(c)
InkaStats Data Science Solutions | Medical Branch
@github/psotob91
A menudo tenemos interés en describir los datos relacionados al tiempo de ocurrencia de algún evento de interés.
Los eventos de interés pueden ser muchos: muerte, cura, desarrollo de alguna enfermedad, entre abandono al tratamiento, etc.
Un problema que tenemos es que no siempre podemos observar el tiempo de ocurrencia del evento:
Los eventos pueden tardar mucho tiempo y el estudio acaba antes de que todos tengan el evento.
Los participantes abandonan el estudio.
Idealmente, todos los participantes desarrollan el evento de interés.
Podemos estimar entonces:
\[\bar{x}_{tiempo} = \frac{10+8+10+6+10}{5} = 8.8 \text{ meses}\]
\[IA(t = 6 \text{ meses}) = 100 \times \frac{1}{5} = 20\%\]
\[IA(t = 8 \text{ meses}) = 100 \times \frac{2}{5} = 40\%\]
\[IA(t = 10 \text{ meses}) = 100 \times \frac{5}{5} = 100\%\]
Los tiempos al evento están disponibles para algunos sujetos.
Promediar estos tiempos no es válido
Todavía podemos hallar la IA de manera directa.
Media con los datos observados con tiempos incompletos:
IA de muerte a 6 meses (%):
\[IA(t = 6 \text{ meses}) = 100 \times \frac{1}{5} = 20\%\]
\[IA(t = 8 \text{ meses}) = 100 \times \frac{2}{5} = 40\%\]
IA de muerte a 10 meses (%):
IA no siempre puede calcularse de proporción observada en cada tiempo t.
Hay censuras por pérdida de seguimiento.
No todos los individuos censurados tienen misma oportunidad de ser seguidos en todo el estudio.
No todos los individuos censurados tuvieron mismo tiempo de seguimiento.
Debemos usar algún método tenga en cuenta los tiempos de seguimiento desiguales en los indviduos censurados.
Densidad de incidencia
Tasa de incidencia
Incidencia y Sobrevidad Acumuladas1
Tablas de vida2
Hazards
Nosotros solo nos centraremos en dos medidas que se usan tradicionalmente para describir datos de supervivencia:
Densidad/Tasa de incidencia
Incidencia y Sobrevida Acumulada1
Tabla de vida1
\[\text{Densidad de Incidencia} = \frac{\text{Número de eventos nuevos}}{\text{Total de tiempo-persona}}\]
\[\text{Densidad de Incidencia} = \]
\[\frac{2 \text{ casos}}{9 +8+5+1+14 \text{ meses-persona}}\]
\[= 0.04 \text{ casos por mes-persona}\]
\[= 4 \text{ casos por 100 meses-persona}\]
Medida de la frecuencia con la que un evento de interés ocurre a lo largo de un periodo de tiempo especificado.
Algunos autores diferencian el concepto de densidad del de tasa.
Skzlo, et al., en su libro “Epidemiology beyond the basics”, hace tal distinción e indica lo siguiente:
La densidad se calcula con los datos de individuales de los tiempos de seguimiento de cada persona (tiempo en riesgo).
La tasa se calculan cuando no se cuentan con datos de los tiempos individuales; en cambio, se disponen de datos agregados.
Densidad de incidencia y tasa de incidencia tienen las misma interpretación y supuestos.
Si tenemos datos con los tiempos de seguimiento, es relativamente fácil calcular.
Veamos los datos:
sex age crea eGFR_ckdepi grf_cat time5y eventd5y
1 Female 72 1.60 31.67301 G3b 3.359343 0
2 Female 84 1.98 22.49146 G4 5.000000 0
3 Female 67 1.10 51.55401 G3a 5.000000 0
4 Male 75 1.36 50.47194 G3a 3.589322 0
5 Female 71 1.24 43.41346 G3b 3.764545 0
6 Female 74 1.06 51.60754 G3a 3.657769 0
eventd5ylab death5y
1 Alive w/o Kidney Failure 0
2 Alive w/o Kidney Failure 0
3 Alive w/o Kidney Failure 0
4 Alive w/o Kidney Failure 0
5 Alive w/o Kidney Failure 0
6 Alive w/o Kidney Failure 0
Variables de interés para cálculo de densidad de incidencia:
Cálculo de densidad de incidencia:
0.0492 caso por 1 año-persona
4.92 casos por 100 años-persona
Muchas veces, los datos individuales no tienen los tiempos de seguimiento, sino las fechas.
# A tibble: 6 × 4
case_id date_onset date_outcome outcome
<chr> <chr> <chr> <chr>
1 5fe599 2014-05-13 <NA> <NA>
2 8689b7 2014-05-13 2014-05-18 Recover
3 11f8ea 2014-05-16 2014-05-30 Recover
4 b8812a 2014-05-18 <NA> <NA>
5 893f25 2014-05-21 2014-05-29 Recover
6 be99c8 2014-05-22 2014-05-24 Recover
character a date:library(lubridate)
ebola_data_small <- ebola_data_small %>%
mutate(date_onset = ymd(date_onset),
date_outcome = ymd(date_outcome))
head(ebola_data_small)# A tibble: 6 × 4
case_id date_onset date_outcome outcome
<chr> <date> <date> <chr>
1 5fe599 2014-05-13 NA <NA>
2 8689b7 2014-05-13 2014-05-18 Recover
3 11f8ea 2014-05-16 2014-05-30 Recover
4 b8812a 2014-05-18 NA <NA>
5 893f25 2014-05-21 2014-05-29 Recover
6 be99c8 2014-05-22 2014-05-24 Recover
ebola_data_small <-
ebola_data_small %>%
mutate(tseg_dias = as.double(date_outcome - date_onset))
head(ebola_data_small)# A tibble: 6 × 5
case_id date_onset date_outcome outcome tseg_dias
<chr> <date> <date> <chr> <dbl>
1 5fe599 2014-05-13 NA <NA> NA
2 8689b7 2014-05-13 2014-05-18 Recover 5
3 11f8ea 2014-05-16 2014-05-30 Recover 14
4 b8812a 2014-05-18 NA <NA> NA
5 893f25 2014-05-21 2014-05-29 Recover 8
6 be99c8 2014-05-22 2014-05-24 Recover 2
ebola_data_small <-
ebola_data_small %>%
mutate(outcome2 = case_when(outcome == "Recover" ~ 0,
outcome == "Death" ~ 1,
TRUE ~ as.numeric(NA)))
head(ebola_data_small)# A tibble: 6 × 6
case_id date_onset date_outcome outcome tseg_dias outcome2
<chr> <date> <date> <chr> <dbl> <dbl>
1 5fe599 2014-05-13 NA <NA> NA NA
2 8689b7 2014-05-13 2014-05-18 Recover 5 0
3 11f8ea 2014-05-16 2014-05-30 Recover 14 0
4 b8812a 2014-05-18 NA <NA> NA NA
5 893f25 2014-05-21 2014-05-29 Recover 8 0
6 be99c8 2014-05-22 2014-05-24 Recover 2 0
Recordar que hay NA
Si hay NA, siempre hay que considerarlo
ebola_data_small %>%
summarise(DI = sum(outcome2, na.rm = TRUE) / sum(tseg_dias, na.rm = TRUE) * 100)# A tibble: 1 × 1
DI
<dbl>
1 4.75
ncasos region pop_mid_year
1 1245 Lima 214555
2 5234 Arequipa 341421
3 2312 Ayacucho 54545
4 4524 Tacna 123454
tasa_incidencia
1 18.14094
La tasa de incidencia de la enfermedad para todas las regiones estudiadas fue de 18.14 casos por 1000 años-persona.
ncasos region pop_mid_year tasa_incidencia
1 1245 Lima 214555 5.802708
2 5234 Arequipa 341421 15.330047
3 2312 Ayacucho 54545 42.387020
4 4524 Tacna 123454 36.645228
La tasa de incidencia de la enfermedad para Lima fue de 5.8 casos por 1000 años-persona. En Arequipa fue de 15.3 casos por 1000 años-persona. La TI de la enfermedad en Ayacucho fue de 42.4 casos por 1000 años-persona y en Tacna fue de 36.6 casos por 1000 años-persona.
La incidencia acumulada en un tiempo t o \(IA(t)\) o proporción de incidencia en un tiempo t, es la probabilidad de desarrollar un evento nuevo en un intervalo t de tiempo.
\[IA(t) = \frac{\text{Nº de eventos nuevos hasta tiempo t}}{\text{Total de individuos en riesgo al incio del estudio}}\]
Si hay pérdidas de seguimiento, debe usarse algún método que tenga en cuenta estas censuras:
El método de Kaplan y Meier es una opción popular en Bioestadística y Epidemiología.
Otros métodos también existen, por ejemplo, el estimador de Nelson-Aalen.
La sobrevida acumulada en tiempo t o \(SA(t)\) es la probabilidad de que un individuo sobreviva más allá de un tiempo \(t\).
\[SA(t) = \frac{\text{Nº de individuos aún en riesgo en tiempo t}}{\text{Total de individuos en riesgo al incio del estudio}}\]
Con pérdida de seguimiento, se puede usar el método de KM.
La \(SA(t)\) es el complemento de la \(IA(t)\):
\[SA(t) = 1 - IA(t)\]
También conocida como método de producto-límite, permite estimar la \(SA(t)\) y la \(IA(t)\) teniendo en cuenta que hay censuras.
No entraremos en detalles técnicos.
KM permite estimar cada \(SA(t)\) e \(IA(t)\) en los puntos de tiempo \(t\) obervados.
Usaremos los datos de ébola:
No tenemos una sola variable de interés observada, tenemos un vector de dos variables:
Tiempo de seguimiento observado:
Indicador de evento y censura:
Importante que el indicador de evento/censura sea numérico con 1 y 0.
Call: survfit(formula = Surv(tseg_dias, outcome2) ~ 1, data = ebola_data_small)
2389 observations deleted due to missingness
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 3499 30 0.991 0.00156 0.988 0.994
2 3464 69 0.972 0.00281 0.966 0.977
3 3385 149 0.929 0.00435 0.920 0.937
4 3216 194 0.873 0.00565 0.862 0.884
5 2989 214 0.810 0.00667 0.797 0.824
6 2747 210 0.748 0.00740 0.734 0.763
7 2470 179 0.694 0.00790 0.679 0.710
8 2232 167 0.642 0.00827 0.626 0.659
9 1992 145 0.595 0.00853 0.579 0.612
10 1772 109 0.559 0.00870 0.542 0.576
11 1592 119 0.517 0.00885 0.500 0.535
12 1398 89 0.484 0.00895 0.467 0.502
13 1214 55 0.462 0.00902 0.445 0.480
14 1097 43 0.444 0.00908 0.427 0.462
15 989 31 0.430 0.00913 0.413 0.448
16 895 48 0.407 0.00923 0.389 0.426
17 775 29 0.392 0.00931 0.374 0.411
18 698 21 0.380 0.00938 0.362 0.399
19 616 7 0.376 0.00941 0.358 0.395
20 554 4 0.373 0.00944 0.355 0.392
21 497 13 0.363 0.00957 0.345 0.383
22 427 8 0.357 0.00969 0.338 0.376
23 378 5 0.352 0.00979 0.333 0.372
24 346 4 0.348 0.00989 0.329 0.368
25 305 4 0.343 0.01002 0.324 0.363
26 274 3 0.339 0.01014 0.320 0.360
27 242 1 0.338 0.01019 0.319 0.359
29 190 1 0.336 0.01029 0.317 0.357
38 60 1 0.331 0.01155 0.309 0.354
El problema con el paquete survival() es que usa R base, por lo que es difícil de manipualr.
El paquete ggsurvfit() utiliza internamente survival() pero permite trabajar con R tidy potenciando el reporte de tablas de vida.
Supervivencia acumulada
surv_ebola <- survfit2(Surv(tseg_dias, outcome2) ~ 1, data = ebola_data_small)
tidy_survfit(surv_ebola)# A tibble: 59 × 14
time n.risk n.event n.censor cum.event cum.censor estimate std.error
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 3499 0 0 0 0 1 0
2 1 3499 30 5 30 5 0.991 0.00157
3 2 3464 69 10 99 15 0.972 0.00289
4 3 3385 149 20 248 35 0.929 0.00468
5 4 3216 194 33 442 68 0.873 0.00647
6 5 2989 214 28 656 96 0.810 0.00823
7 6 2747 210 67 866 163 0.748 0.00989
8 7 2470 179 59 1045 222 0.694 0.0114
9 8 2232 167 73 1212 295 0.642 0.0129
10 9 1992 145 75 1357 370 0.595 0.0143
# … with 49 more rows, and 6 more variables: conf.high <dbl>, conf.low <dbl>,
# estimate_type <chr>, estimate_type_label <chr>, monotonicity_type <chr>,
# conf.level <dbl>
Incidencia acumulada
# A tibble: 59 × 14
time n.risk n.event n.censor cum.event cum.censor estimate std.error
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 3499 0 0 0 0 0 0
2 1 3499 30 5 30 5 0.00857 0.00157
3 2 3464 69 10 99 15 0.0283 0.00289
4 3 3385 149 20 248 35 0.0711 0.00468
5 4 3216 194 33 442 68 0.127 0.00647
6 5 2989 214 28 656 96 0.190 0.00823
7 6 2747 210 67 866 163 0.252 0.00989
8 7 2470 179 59 1045 222 0.306 0.0114
9 8 2232 167 73 1212 295 0.358 0.0129
10 9 1992 145 75 1357 370 0.405 0.0143
# … with 49 more rows, and 6 more variables: conf.high <dbl>, conf.low <dbl>,
# estimate_type <chr>, estimate_type_label <chr>, monotonicity_type <chr>,
# conf.level <dbl>
Tiempos específicos de interés
# A tibble: 3 × 14
time n.risk n.event n.censor cum.event cum.censor estimate std.error
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 3499 30 5 30 5 0.00857 0.00157
2 2 3464 69 10 99 15 0.0283 0.00289
3 5 2989 557 81 656 96 0.190 0.00823
# … with 6 more variables: conf.high <dbl>, conf.low <dbl>,
# estimate_type <chr>, estimate_type_label <chr>, monotonicity_type <chr>,
# conf.level <dbl>
Selecciona solo columnas de interés
Más información sobre análisis de datos de supervivencia:
{survival}: https://cran.r-project.org/web/packages/survival/index.html (Ver Sección Vignettes)
{ggsurvfit}: https://www.danieldsjoberg.com/ggsurvfit/index.html
Un tutorial básico para epidemiólogos sobre análisis de supevivencia: https://epirhandbook.com/en/survival-analysis.html
Una lista amplia de paquetes en R para análisis de supervivencia: https://cran.r-project.org/web/views/Survival.html
Un (más avanzado) tutorial interesante sobre análisis de supervivencia: https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html
Es la versión gráfica de la tabla de vida.
Hay varios métodos para obtenerla.
Veremos algunas funciones para obtener gráficos para reporte.
Solo nos centraremos en el método no paramétrico de Kaplan-Meier para describir datos de supervivencia.
Call: survfit(formula = Surv(tseg_dias, outcome2) ~ 1, data = ebola_data_small)
2389 observations deleted due to missingness
time n.risk n.event survival std.err lower 95% CI upper 95% CI
1 3499 30 0.991 0.00156 0.988 0.994
2 3464 69 0.972 0.00281 0.966 0.977
3 3385 149 0.929 0.00435 0.920 0.937
4 3216 194 0.873 0.00565 0.862 0.884
5 2989 214 0.810 0.00667 0.797 0.824
6 2747 210 0.748 0.00740 0.734 0.763
7 2470 179 0.694 0.00790 0.679 0.710
8 2232 167 0.642 0.00827 0.626 0.659
9 1992 145 0.595 0.00853 0.579 0.612
10 1772 109 0.559 0.00870 0.542 0.576
11 1592 119 0.517 0.00885 0.500 0.535
12 1398 89 0.484 0.00895 0.467 0.502
13 1214 55 0.462 0.00902 0.445 0.480
14 1097 43 0.444 0.00908 0.427 0.462
15 989 31 0.430 0.00913 0.413 0.448
16 895 48 0.407 0.00923 0.389 0.426
17 775 29 0.392 0.00931 0.374 0.411
18 698 21 0.380 0.00938 0.362 0.399
19 616 7 0.376 0.00941 0.358 0.395
20 554 4 0.373 0.00944 0.355 0.392
21 497 13 0.363 0.00957 0.345 0.383
22 427 8 0.357 0.00969 0.338 0.376
23 378 5 0.352 0.00979 0.333 0.372
24 346 4 0.348 0.00989 0.329 0.368
25 305 4 0.343 0.01002 0.324 0.363
26 274 3 0.339 0.01014 0.320 0.360
27 242 1 0.338 0.01019 0.319 0.359
29 190 1 0.336 0.01029 0.317 0.357
38 60 1 0.331 0.01155 0.309 0.354
Podemos usar la función ggsurvfit() para realizar gráficos publicables de la KM.
El paquete {survminer} también ofrece funciones para gráficos ggplot2().
Sin embargo, los objetos que crean también tienen atributos adicionales:
Personalmente creo que esto genera algunas inconsistencias y obliga a salirse un poco de la sintaxis {ggplot2}
No veremos {survminer} en detalle, porque ggsurvfit() lo supera por creces.
Solo veremos un par de ejemplos para conocer, por cultura general, el paquete en general:
Más información sobre gráficos KM para datos de supervivencia:
{ggsurvfit}: https://www.danieldsjoberg.com/ggsurvfit/index.html
Gráficos: https://www.danieldsjoberg.com/ggsurvfit/articles/gallery.html
Temas: https://www.danieldsjoberg.com/ggsurvfit/articles/themes.html
Listado de todas las funciones: https://www.danieldsjoberg.com/ggsurvfit/reference/index.html
{survminer}: https://rpkgs.datanovia.com/survminer/index.html
Un tutorial para epidemiólogos sober survminer: https://epirhandbook.com/en/survival-analysis.html
Abra el proyecto var_tiempo.Rproj y dentro de este, abra el archivo quarto var_tiempo.qmd.
Siga las instrucciones indicadas en este.
Renderice el archivo quarto final.
10:00
Datos de tiempo de supervivencia